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Abstract. This paper presents estimates of the convergence rate and complexity of an algebraic 
multilevel preconditioner based on piecewise constant coarse vector spaces applied to the graph 
Laplacian. A bound is derived on the energy norm of the projection operator onto any piecewise 
constant vector space, which results in an estimate of the two-level convergence rate where the 
coarse level graph is obtained by matching. The two-level convergence of the method is then used 
to establish the convergence of an Algebraic Multilevel Iteration that uses the two-level scheme 
recursively. On structured grids, the method is proven to have convergence rate ~ (1 — 1/logn) 
and 0(n log n) complexity for each cycle, where n denotes the number of unknowns in the given 
problem. Numerical results of the algorithm applied to various graph Laplacians are reported. It 
is also shown that all the theoretical estimates derived for matching can be generalized to the case 
of aggregates containing more than two vertices. 



1. Introduction 

Algebraic Multigrid (AMG) attempts to mimic the main components of Geometric Multigrid in 
an algebraic fashion, that is, by using information from the coefficient matrix only to construct the 
multilevel solver. The basic algorithm uses a setup phase to construct a nested sequence of coarse 
spaces that are then used in the solve phase to compute the solution. The two main approaches to 
the AMG setup algorithm are classical AMG [3 i] and (smoothed) aggregation AMG [ISj IT71I91 
[m [11] , which are distinguished by the type of coarse variables used in the construction of AMG 
interpolation. 

In the classical AMG algorithm, the coarse variables are chosen using a coloring algorithm 
which is designed to find a suitable maximal independent subset of the fine variables. Then, given 
the coarse degrees of freedom, a row of interpolation is constructed for each fine point from its 
neighboring coarse points. In contrast, the aggregation-based AMG setup algorithm partitions the 
fine variables into disjoint subdomains, called aggregates. Then, a column (or several columns as 
in |18j ) of interpolation is associated to each aggregate, which has nonzero entries only for the 
unknowns belonging to this aggregate. The focus of this paper is on the development and, in 
particular, the analysis of the latter aggregation-type methods. 

The idea of aggregating unknowns to coarsen a system of discretized partial differential equations 
dates back to work by Leont'ev in 1959 [12]. Simon and Ando developed a related technique for 
aggregating dynamic systems in 1961 [13j and a two- grid aggregation-based scheme was consid- 
ered in the context of solving Markov chain systems by Takahashi in 1975 [TS]. Aggregation-based 
methods have been studied extensively since and numerous algorithms and theoretical results have 
followed [Hi [T8l [TT] . Vanek introduced an extension of these methods known as smoothed aggre- 
gation multigrid in which smoothing steps are applied to the columns of the aggregation-based 
interpolation operator to accelerate two-level convergence and a modification of this two-level algo- 
rithm with overcorrection is presented in [T7] . A multilevel smoothed aggregation algorithm and its 
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convergence analysis are found in |16j and, in [19j, an improved convergence theory of the method 
is presented. The latter theory is then extended to allow for aggressive coarsening, provided an 
appropriate polynomial smoother is used |8j . A further generalization known as adaptive smoothed 
aggregation is developed in [J]. Variants of the above approaches continue to be developed for 
use in scientific computing and have been developed for higher order partial differential equations 
|18j , convection diffusion problems , Markov chains [1^ [2] , and the Dirac equation in quantum 
chromodynamics [5]. 

In this paper, an aggregation based Algebraic Multigrid method for the graph Laplacian is 
presented. The approach constructs the sequence of coarse graphs recursively using a pair-wise 
aggregation, or matching, form of interpolation. However, it is demonstrated here, that the con- 
vergence rate of a two-level method based on such a construction is uniformly bounded for the 
graph Laplacian on general graphs and, thus, can be used within an Algebraic Multilevel Iteration 
(AMLI) [U [19] as a preconditioner to the Conjugate Gradient iteration to obtain a nearly optimal 
solver. A noteworthy feature of the approach is its simplicity, which makes it possible to analyze 
the convergence and complexity of the method with few assumptions and without any geometric 
information. 

The remainder of the paper is organized as follows. In Section[2| we introduce the graph Laplacian 
problem and discuss some of its applications. In Section [3| we introduce a graph matching algorithm 
and demonstrate that the energy norm of the £2 projection onto the coarse space is a key quantity 
in deriving convergence and complexity estimates of the method. Additionally, we introduce an 
approach computing an approximation of the energy norm of this projection operator. In Section 
|4j we present an analysis on the two-level method for the graph Laplacian operator. In Section [5j 
we consider the convergence and complexity of the resulting AMLI method, and in Section [6] we 
provide numerical results and address some practical issues of the method. 

2. Problem formulation and notation 

Consider an unweighted connected graph Q = (V,<S), where V denotes the set of vertices and 
£ denotes the set of edges of Q. The variational problem considered here is as follows: Given an 
/ satisfying = 0, where 1 is a constant vector, find a u G M", where n = |V| denotes the 

cardinality of the set of vertices V, such that 

(2.1) {Au,v) = {f,v), V-rGR", 
where 

(2.2) {Au,v)= J2 {u,-uj){v,-Vj), (/,^) = ^/^^i , (/,1) = 5^/.. 

k={i,j)e£ iev iev 

Define the discrete gradient operator B : RI^I i— )• RI^I such that 

{Bu)k = Ui-Uj, k = {i,j)££, i<j, 
where ej and ej are standard Euclidean bases. The operator A is named graph Laplacian since 

A = B'^B. 

The operator A is symmetric and positive semi-definite and its kernel is the space spanned by 
the constant vector. These properties can also be verified by the matrix form of A defined in the 
following way 



di i = j; 

i / j,(«,j) ^ 
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where di is the degree of the i-th vertex. 

Efficient multilevel graph Laplacian solvers are important in numerous application areas, in- 
cluding finite element and finite difference discretizations of elliptic partial differential equations 
(PDE), data mining, clustering in images, and as a preconditioner for weighted graph Laplacians. 
Moreover, the theory developed here for multilevel aggregation solvers applied to graph Laplacians 
should provide insights on how to design a solver for more general weighted graph Laplacians, which 
cover also anisotropic diffusion problems. Two generalizations of the graph Laplacian systems are 
as follows. 

• Weighted graph Laplacians: Assume that the graph is weighted and the A;-th edge is 
assigned a weight Wk , then the corresponding bilinear form of A is 

{Au, v) = ^ Wk{ui - Uj){vi - Vj). 

k={i,j)e£ 

Define D : RI^I ^ RI^I as a diagonal matrix whose fe-th diagonal entry is equal to Wk, then 
the matrix A can be decomposed as 

A = B^DB. 

Finite element and finite difference discretizations of elliptic PDEs with Neumann boundary 

conditions results in such weighted graph Laplacians. 

• Positive definite matrices: Assume that A is defined in terms of the bilinear form 

{Au, v)= ^ {ui - Uj){vi - Vj) + ^ UiVi = (AgU, v) + {Atu, v). 

k={i,j)e£ iev 

By introducing a Lagrange multiplier y, the system Au = f can be rewritten as an aug- 
mented linear system 

(A, + At -Atl\ (u\^( f \ 

[-l^A, l^AlJ [yj [-l^fj' 

These graph Laplacians with lower order terms are similar to discretized PDE with Dirichlet 
boundary conditions. A solution for this augmented linear system directly results to the 
solution of Au = f. 

The present paper focuses on designing a multilevel preconditioner that is constructed by apply- 
ing recursively a space decomposition based on graph matching. The aim is to analyze the matching 
AMLI solver for the graph Laplacian in detail as a first step in gaining an in-depth understanding 
of a multilevel solver elliptic PDEs. The extension of the proposed algorithm to general graph 
problems is also a subject of current research. 

3. Space decomposition based on matching 

In this section, an outline of the basic idea of matching is provided, a commutative diagram 
which can be used to estimate the energy norm of the £2 projection onto the piece- wise constant 
coarse vector space resulting from a matching in a graph is given, and some auxiliary results that 
are needed later on in the convergence analysis are discussed. 

3.1. Subspaces by graph partitioning and graph matching. A graph partitioning of ^ = 
{V,6) is a set of subgraphs Qi = {Vi,6j) such that 

In this paper, all subgraphs are assumed to be non empty and connected. The simplest non trivial 
example of such a graph partitioning is a matching, i.e, a collection (subset A^) of edges in £ such 
that no two edges in Ai are incident. 
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For a given graph partitioning, subspaces of V = M'^I are defined as 

Vc = {v G V\ V = constant on each Vi }. 

Note that each vertex in Q corresponds to a connected subgraph S ofQ and every vertex of Q belongs 
to exactly one such component. The vectors from Vc are constants on these connected subgraphs. 
Of importance is the £2 orthogonal projection on Vc, which is denoted by Q, and defined as follows: 

(3.1) (Qv), = -L ^ G Vfe. 

Given a graph partitioning, the coarse graph Qc = { Vc, £c} is defined by assuming that all vertices 
in a subgraph form an equivalence class, and that Vc and £c are the quotient set of V and £ under 
this equivalence relation. That is, any vertex in Vc corresponds to a subgraph in the partitioning 
of G, and the edge exists in £c if and only if the i-th and j-th subgraphs are connected in the 
graph Q. Figure [T] is an example of matching of a graph and the resulting coarse graph. 




Figure 1. Matching on a graph Q (left) and the coarse graph Qc (right) 



As mentioned above, the reason to focus on matching is that it simplifies the computation of 
several key quantities used in the upcoming estimates derived for a perfect matching and it is 
possible to show that a matching which is not perfect can be analyzed in a similar way. 

3.2. Commutative diagram. Let B be the discrete gradient of a graph Laplacian A, as defined 



in (2.1), and Q be defined as in (3.1). Assume that there exists an operator 11^ such that the 
following commutative diagram holds true: 

m|V| ^\£\ 
n 

Vc > mI^i 

B 

The proof of this assumption is provided later on. From the commutative relation BQ = IIB it 
follows that 

(3.2) \Qv\\ = \\BQvf = \\UBvf < ||n|| VIa- 

Thus, an estimate on the A-semi-norm of Q amounts to an estimate of the ^2 norm of H. In the 
next subsection, an explicit form of n is constructed and an estimate of its £2 norm is derived. 

Remark 3.1. A more general approach for weighted graph Laplacians is to assume that the weight 
matrix D ^ I, therefore the bound on the norm \Q\a becomes 

\Qv\\ = {DBQv,BQv) = {DUBv,UBv) < \\D^/'^nkD-'^/'^f\v\\, 

where D can have some negative weights, which results in a matrix D^^^H^.D^^^'^ this is complex 
valued. A detailed analysis in such a setting and the application of this idea to anisotropic diffusion 
problems are discussed in 
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3.3. Construction of 11 in case of piece-wise constant spaces. Here, we proceed with an 
explicit construction and £2 norm estimate of tiie operator 11. 

For any grapli partitioning in which the subgraphs are connected, a given edge belongs to the set 
of "internal edges" , whose vertices belong to the same subgraph, or to the set of "external edges" , 
whose vertices belong to two distinct subgraphs. For example, let Gi and G2 denote the subgraphs 
1 and 2 in Fig. [2| then ki is an internal edge and k2 is an external edge. 



Figure 2. Connected components and the construction of 




Since the vector Qv has the same value on the two endpoints of the edge ki, we have that 
{BQv)k^ = 0. Accordingly, all entries in (n)^^, the /ci-th row of 11, are set to zero: 

iUBv)k, = iU)k,Bv = 0. 

For the external edge k2, it follows that (n)^^ satisfies 

(3.3) (n)fe(S^) = {BQv)k, = ^ E - E ^^2, 

for every v. The following Lemma is useful in computing explicitly the entries of (11)^2 • 

Lemma 3.2. Let A : M" 1— t- M" be a positive semidefinite operator and let {xi}f=i be a basis in M"". 
Assume that the null space of A is one dimensional, namely there exist a nonzero vector s such 
that Ker(^) = span(s), and for every integer 1 < i < n we have {xi, s) = 1. We then have: 

(i) For any i, the operator A : M" 1— t- M" with Au = {Au + (xi, u)xi) is invertible. 

(ii) The following identity holds for all u G M"; 

1 1 ~_i 

-{u,s) - {u,Xi) = AA s,Au). 



(s,s) ' ' (s,s) 



Proof. To establish (i) it suffices to show that Av = implies u = 0. Assuming that Av = for 

= {Av, v) = {Av, v) + {xi,vf. 

Note that both terms on the right side of the above identity are nonnegative and, hence, their sum 
can be zero if and only if both terms are zero. Since A is positive semidefinite by assumption with 
one dimensional null space, from {Av,v) = we conclude that v = as for some a G M. For the 
second term we have that = {Xi^v)"^ = Q;^(xi,'S)^, and since {xi^s) / for all i, it follows that 



a = and hence v = Q. This proves (i) 
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Now, applying I (i) I the result [(li)] follows: 



{s,s) ' {s,s) ' ' (s,s) 

^ -{A-h, Au) - j^ixi, u){s, A-^x^) 



{s,s) ' (s,s) 

{s,u) - r{xi,u){s,A-^Xi) 



{s,s) ' (s,s) 
{s,s) 

Here, the equality As = As + {xi, s)xi = Xi was used, implying that A~^Xi = s . □ 

Remark 3.3. A special case is given by taking s = 1 and Xi = which denote the standard 
Euclidean bases. Then, it follows that 

(3.4) Ui = {u)--{{A + eief)-^l,Au) 



n 



in which {u) := - Yl^=i '^i denotes the average value of 



u. 



Next, denote by Bm the restriction of B to a subgraph Qmi and set Am '■= BmBm- Then, ui can 
be expressed as the l-ih component of u for u in (3.5). 

(3.5) Ui = {u)m + |T^(-Bm(^m + e;ef )~4m, BmlmU), 

where 1 = 1... \ Vm\ are the local indices of the vertex set Vm, the operator {■)m and the term Im 
are the averaging operator and the constant vector restricted on the subgraph Qm-, and : MI^I i— )• 
rI^™! maps the global edge indices to the local edge indices. 

Applying this formula for Q\ and Q2, gives the row of the operator 11 on the edge k2 that connects 
Qi and Q2 as follows 

(3.6) (n),, = CjTj + el - CM- 
Here Ci is given by 

I l^il 



which then makes the summation in ( |3.6| valid. The vector C2 is defined in a similar way. 

Assume that the global indices of the vertices in Qi and Q2 are ordered consecutively as decreasing 
integers starting at A;2 — 1 and increasing integers starting at ^2 + 1. Then, the /c2-th row of H can 
be expressed as 

(3.7) (n)fc, = [0, . . . , 0, Cf, 1, C^, 0, . . . , 0] 



where the number 1 is on the A;2-th position in this row of IT. Note that from (3.5) it follows that 
the property (3.3) holds for this construction of H, since by definition of Q 

{BQv)k, = {v)i-{v)2 

= '"i- + + eief)~^li,BiVi) - T^{B2{A2 + ejej)-^l2, B2V2) 

= {ek„Bv) + {Ci,Biv) + {C2,B2v) 

where k2 = and i and j, both in local indices, are the incident vertices of /c2. 
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4. A TWO-LEVEL METHOD 



In this section, the ^2-orthogonal projection given in (3.1 ) based on a matching A4 is proven to be 
stable assuming the maximum degree of the graph Q is bounded. Then, a two-level preconditioner is 
derived and the condition number of the system preconditioned by this two-level method is proven 
to be uniformly bounded (under the same assumption). 

4.1. Two- level stability. The construction of 11 for a matching Ai proceeds as follows. First, 
note that all rows of 11 that correspond to an edge k = G M are identically zero. On the 

other hand, if the edge k = ^ M, then it is an external edge and, thus, by (3.7), the k-th. row 
of n is 



(n)* 



[0,... 




,0 



Hence, 
(4.1) 



(n) 



kl 







k ^ M and / = k; 

k ^ M,l e M and Ink ^ 

elsewhere. 



The alternative way of describing the entries in 11 is by showing that, 

' 1 I ^ M and k = I; 



(4.2) 



(n) 



kl 







I £ M,k ^ M and knl 
elsewhere. 



Formula (4.1) implies that, the A;-th row of 11 can be a zero row if /c G A4, or a row with 3 

1 + I ± 1/2| + I ± 1/2| 



non-zero entries if k ^ Ai, which results to 



n 



max 

k 



J2\^ki\ 



2. 



Formula (4.2 ) implies that, the l-th. column of IT can have exactly 1 non-zero entry if / ^ M, or s 
non-zeros entries whose values are ±1/2 if / G M. Here s is the number of edges satisfying k ^ M 
and A; n / 7^ for any given I G M, thus is bounded by 2d — 2, where d is the maximum degree of 
the graph, since an edges can have at most 2d — 2 neighboring edges. This leads to 

||n||i = max ^iHfczl =max(l,(2d-2)| ±1/2|) =max(l,d-l). 

' k 

On a graph whose maximal degree is larger or equal to 2, the estimates on the infinity norm and 
£i norm of H result to the following estimate on /9(Hn^): 



(4.3) 



p(ntf 



ml < 



nililini 



2d 



Remark 4.1. Applying Gerschgorin's theorem directly to the matrix leads to a sharper esti- 

mate: /9(Hn'^) < d. 



Formula (|4.3|) implies directly the following lemma. 
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Lemma 4.2. On any graph whose maximum degree is 2 (e.g. such graph is a path), the operator 
n defined in (4.1) satisfies Ilk{Bv) = [BQv)k and the following estimate holds 

IQIi < ||n||^ < 2d - 2 = 2. 

Numerical tests show that this is a sharp estimate on the semi- norm \Q\a and that leads to fast 
convergent and reliable AMG methods. 

4.2. A two-level preconditioner. Here, using an estimate of the stability of the matching projec- 
tion (i.e, the norm \Q\a^ where Q is defined via the matching) two- level convergence is established. 
Assume that for a graph Laplacian A : i— )• M" a perfect matching is given and consider the 
n X n/2 matrix P whose fc-th column is given by 

(4.4) (P)fc = ei, +e,,, 

where A; = 1, ...,n/2 and {ik,jk) is the k-th edge in M. Further, define Q to be the £2 projection 
from M" to {Pv\v G M"/^}, i.e., 

Q = P{P'^P)-'^P^. 

Similar to the definition of P, define Y as the n x n/2 matrix whose columns are given by 

(4.5) iY)k = ei^-ej^, 

where k = l,...,n/2 and {ik,jk) is the k-th edge in A4. Then, the matrix (Y,P) is orthogonal and 
the columns of Y and P form a hierarchical bases, which can be used to relate the two-level method 
to a block factorization as follows. 
Given A, P, and Y, define 



Y'^AY Y'^AP^ 
P^AY P^AP. 




A = {Y,PfA{Y,P) 
A direct calculation then shows that 

A = L 
where 

(4.6) S = P'^AP - P'^ AY (Y'^ AY)-^Y'^ AP 

is the Schur complement and 

P^AYiY^AY)-^ I 



(4.7) L 



Next, define Qc as the unweighted coarse graph and denote by Ac the graph Laplacian of Qc- In 
contrast to most of the existing AMG methods, here Ac P^AP, except in special cases, e.g., for 
1 dimensional problems. Let, cr be a positive constant such that 

/ , (APv,Pv) 
(4.8) a= sup \, / . 

v:iv,l)=0 {AcV,v) 

Then, the fact that all weights in the graph corresponding to P^AP are larger than or equal to 
one implies {APv,Pv) > {AcV,v),\/v, and 

(aAcV,v) , , , , 

^ ' ^ G [I, a], yv : {v,l) = 0. 



{APv,Pv) 
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Consider the two-level preconditioner G which uses the coarse graph Laplacian Ac by 

\ aAc) 

Let M be a preconditioner for Y^AY, and D be a preconditioner for the graph Laplacian Ac- 
Then, a two-level preconditioner B is defined by 

(4.9) g^-/'A/(M + M--r-AF)-'M- o\- 

\^ aDj 

where 

L-( ' ") 

yp^AYM-^ I J 

As observed in |10j and [19], this gives a block matrix representation of the two-level method 

/ - {Y, P)G\Y, PfA = (I - Y{Y^AY)-^Y^A){I - P{aAc)^ P'^ A){I - Y{Y^ AY)-^Y^ A) 
I - {Y,P)B\Y,PfA = {I -YM''^Y'^A){I - P{aD)^P^A){I -YM-^Y'^A), 

where the pseudo-inverse operator denoted by is used since the graph Laplacian is semi-definite. 
The convergence of the two-level method can now be estimated by comparing A and the precondi- 
tioner B. 

The remainder of this section is dedicated to establishing a spectral equivalence between A and 
B for the two-level matching algorithm. The proof uses the following Lemma. 



Lemma 4.3. For any x G iR"/^ the Schur complement S as given in {4-6) satisfies 
(4.10) {Sx, x) = inf (A{Yw + Px), {Yw + Pa;)) . 

w 

Proof. Note that 

{AY(Y^AY)-^Y^APx,Px) = {AY{Y^ AY)-^Y'^ APx,Y{Y'^ AY)-^Y^ APx) 

= \\Y{Y^AY)-^Y^APx\\\, 

because here, Y {Y^ AY)~^Y^ APx is the A orthogonal projection of Px onto the space spanned 
by the columns of Y and, thus, minimizes the distance (in A norm) between Px and this space. 
Hence, 

iSx,x) = \\Px\\\ - \\Y{Y^AY)-^Y^APx\\\ 

= inf (A(Yw + Px),(Yw + Px)) □ 

Let 1 be a vector satisfying (y, P)\ = 1, then the following lemma now holds. 



Lemma 4.4. Let Cg = a\QW, where a is defined as in (4.8). Then for any v, such that (v, 1) = 0, 
we have 

(4.11) 7m 

[Av, v) 



Proof. By Lemma |4.3| we have 

{APx, Px) > inf {A{Yw + Px) , {Yw + Px)) . 
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Furthermore, 



{APx, Px) 
{Sx, x) 



{APx,Px) 



inf^ {A{Yw + Px), {Yw + Px)) 
{APx,Px) 



sup 



{A{Yw + Px), {Yw + Px)) 



{AQu,Qu) (AQu,Qu) ^ 
= sup — ^<sup— — = \Q\^. 

Note that the only difference between the preconditioners G and A is that the former matrix uses 
aAc, whereas the latter uses S to define the 2-2 block. The spectral equivalence constant between 
the operators aAc and 5 is obtained as follows: 



inf 



a{AcU,u) .^^ {APv,Pv) ^ a{AcW,w) 

{Sw, w) 



inf 



u {APu,Pu) V {Sv,v) 
which implies 

a{AcW, w) 



< 



< sup 



cr{AcU,u) {APv,Pv) 
{APu, Pu) ""t"^ {Sv,v) ' 
yw : (w,l) = 0, 



(Sw, w) 



e[l,a\Q\l], yw:iw,l) = 0. 



Hence, for any x and y 

T 





{AYx,Yx) + {APy,Py) 
{AYx,Yx) + {Sy,y) 



which is equivalent to (4.11) since L is nonsingular. 



□ 



Since the two- level method G requires exact solvers for Y^ AY and the graph Laplacian Ac, the 
convergence rate of a method that uses B which is defined by replacing these exact solves with 
approximate ones is of interest. Combining Lemma 4.4 and the two-level convergence estimate 
(Theorem 4.2 in [lOj), yields the following result. 

Theorem 4.5. // the preconditioners M and D are spectrally equivalent to Y^ AY and Ac such 
that 

({M^ + M -Y^AY)-^Mu,Mu) ^ ^ , (Dw,w) , , , , 

'-e[l,Ks] and }- — ^G[l,r/], Vw, w : (w, 1) = 0, 



then 
(4.12) 



{AYu,Yu) 

{Bv,v) 



{AcW,w) 



{Av,v) 



G [l,{Ks + av-l)\Q\\], 



\/v : {v,l) = 0. 



Note that this estimate reduces to (4.11) when M = Y AY and D = Ac- 

4.3. Convergence estimate for matching. We here show the sharpness of the estimation given 
by Theorem 4.5 when the graph Laplacian corresponds to a structured grid, and the coarse space 
is given by aligned matching. 

Define an m-dimensional hypercubic grid as the graph Laplacian Q = (V, £) such that the 
following conditions are satisfied. 
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(1) A vertex S V corresponds to an vector v S M™, and {v, ej) G [1,2,..., sj], j = 1, 2, . . . , m. 
Here Cj is an Euclidean basis and si, S2, ■ ■ ■ , Sm are given positive integers that represent 
the numbers of vertices along all dimensions. 

(2) An edge k = {iu,iv) is in the edge set £ if and only u — v = ej and j £ [1, 2, . . . , m]. 
Then the energy norm \Q\a can be estimated for aligned matching on a hypercubic grid A. 

Lemma 4.6. Let Q he an m- dimensional hypercibic grid and k £ [1,2, ... ,m] is a fixed dimension. 
Assume that is an even number. The matching along the k-th dimension is defined as 

M = {I = {iv,iv+ek)\v £ V, and {v,ek) is an odd number }. 

Let Q be the I2 projection onto the piecewise constant space resulting from the matching M . Then 
Q satisfies \Q\a < 2. 

Proof. Define the set il. be the collection of all edges along the k-th dimension, as 

n = {I = {iu,iv)\v - u = Ck}. 

Also define 0, = £\Q and the graph Laplacians ^4^2 and A^, derived from and 0, respectively. □ 

The graphs in the set 0, are paths, whose maximum degree is 2, and M C il. is a also matching 
on these paths. Therefore by Lemma |4.2| it is true that 

(4.13) {AnQu, Qu) < 2{Anu, u). 

On the other hand, the matching is aligned on the set meaning that any two matched pairs are 
connected through or 2 edges in Q, thus the edges in set Q can then be subdivided into many 
sets of edges of the same type, one of which is shown in Fig. |3j Notice that in this figure, the edge 




k t^^^^^W I 

Figure 3. Matching M on a subset of $7 

{i,k) and {j, I) are in il., while {i,j) and {k,l) are in A4. Using the definition of Q, the energy norm 
of Q is estimated on the the subset of indicated by Fig. [s} by 

Ui + Uj Ufc + nA^ 1 2 

' = ^ {{Ui - Uk) + [Uj - Ul)) 
< {Ui - Ukf + {Uj - Ulf . 



Thus implies that 

(4.14) {A^Qu,Qu) <{A^u,u) 



Combining (413) and (4T4) results that {AQu,Qu) < 2{Au,u), or \Q\a < 2. 



Remark 4.7. A similar estimate follows for aligned partitionings consisting of line segments of 
size m. Namely, in this case it can be shown that \Q\\ < m holds. This estimate in turn agrees with 
properties of Chehyshev polynomials, suggesting the use of an AMLL method equipped with certain 
Chebyshev polynomials. Comparing this result with the result from Theorem \4.()\ suggests that using 
a more shape regular partitioning rather than one consisting of lines is more appropriate since this 
gives smaller values of the semi-norm \Q\a- 

A bound on the constant Kg follows by using that AY is well conditioned and that its condition 
number depends on the degree of the graph, but not on the size of the graph. 



12 



J. BRANNICK, Y. CHEN, J. KRAUS, AND L. ZIKATANOV 



Lemma 4.8. Let M he the perfect matching on a graph maximum whose degree is d, and let S he 



defined as in (4.5), then we have 

{AYw,Yw) 



G[4,2d], Vw/O. 



{w, w) 

Proof. The j4-norm of the vector Yw is computed by definition: 

{AYw,Yw)> {iYw)^-{Yw)j)^ = ^ {{Yw)i + {Yw)if = 4w'^ 



k={i,j)eM k={i,j)eM 



We also have 



p{Y^AY) < < ||y^||i||yl||oo||l"||oo = 2d. □ 

From the Lemma it foUows that for any e > there exists a smoother M such that the bound 



on the constant in Theorem 4.5 



is 

Ks <1 + e. 



This result in turn implies that an efficient solver for Y^ AY can be constructed by applying a few 
Conjugate Gradient iterations with an overall cost that is linear with respect to the size of Y^AY. 



The constant a in (4.8) can be estimated by checking the weights of the graph for the graph 
Laplacian P^AP. Taking any two distinct subgraphs (edges) in the matching, say the fc-th and 
l-th. such that A; 7^ /, it follows that the corresponding entry {P^AP)ki is equal to the number of 
exterior edges that connect to these subgraphs. For an aligned matching aligned a fixed dimension 
of a hypercubic grid, these weights are bounded by 2. For any general graph A, the weights in 
P^AP are bounded by 4, since there are at most 4 distinct edges that connect to any other 2 
distinct edges. Then, letting Ac to denote the unweighted graph Laplacian on the graph defined 
by P^AP, and noting that all off-diagonal entries of Ac are equal to —1, it follows that 



a 



for an aligned matching on a hypercubic grid of any dimension; 
for a given matching on any graph. 



Remark 4.9. These estimates can he generalized to other subgraph partitionings in a similar way. 
As an example, consider again a graph for a hypercubic grid of any dimension. Then, for line 
aggregates of size m (aligned with the grid) the following estimate holds 

IQIa — ^) Kg < 1 + e, r] = 1, a < m. 

Such estimates give insight into the design of a nearly optimal multilevel method. Moreover, the 
bounds are sharp enough, namely, the corresponding multilevel method can be proven to have con- 
vergence rate (1 — 1/logn) and 0(n log n) complexity. 



5. Algebraic multilevel iteration (AMLI) based on matching 

In this section, a multilevel method that uses recursively the two-level matching methods from 
Section |4.2| in combination with a polynomial stabilization, also known as Algebraic Multilevel 
Iteration (AMLI) cycle is analyzed. Here, the focus is on proving a nearly optimal convergence 
rate, that is, one which is nearly independent of the number of unknowns n, and at the same time 
has low computational complexity. 
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5.1. Multilevel hierarchy. Assume that Aj = ^ is an n x n graph Laplacian matrix where 
n = 2^ . For k = 1, . . . , J define the matching M.^ and the prolongation operator according to 
(4.4), then compute the graph Laplacian of the coarse graph Qk (Recall that, A^^i ^ P^A^Pk)- 
The index k starts at 1 because the analysis is simpler if the coarsest graph has more than 1 vertex. 
Also, define Y^. and for A^ as in (4.5) and (4.7), and let the two-level preconditioner on each 
level A; be given by 



Gk = Lk{^^ ^y^" ; ]Ll, k = 2,...,J. 



'Y^AkYk 
y aAk-i J 
Then an AMLI preconditioner is defined recursively by 



Bi' = 4, 

V a 'B^_^qk-i{Ak-iBf^_^)J 

= {Yk,PkfB^HYk,Pk), k = 2,...,J, 

where Qkit) is a polynomial on that determines a special coarse level correction on the A:-th level. 
In this case, where an AMLI W-cycle is used, qk{i) is a linear function for all k. 

In the remainder of this section, sufficient conditions for guaranteeing the spectral equivalence 
between the multilevel preconditioner Bj, as defined above, and the graph Laplacian A are derived. 
We first prove two auxiliary results, which are needed in the analysis below. 

Proposition 5.1. Let A : V V and G : V V be symmetric positive semidefinite operators on 
a finite dimensional real Hilhert space V . Suppose that the following spectral equivalence holds: 

(5.1) cq{Av,v) <{Gv,v) <ci{Av,v), cq > 0, ci > 0. 
Then, we also have that 

(5.2) c{^{A^v, v) < {G^v, v) < c^^iA^v, v). 

Proof. Observe that the spectral equivalence given in ( |5.1[ ) implies that A and G have the same 
null-space (and also same range, because they are symmetric). Also, note that, if v is in this null 
space, then ( |5.2[ ) trivially holds. Thus, without loss of generality, we restrict our considerations 
below to V from the range of G and A. 

After change of variables w = (^A^^^^'^v from the upper bound in (5.1) we may conclude that 

-<ci, and hence, \\G^^^ (A^) ^^f < ci. 



Since Gi/2(^t)V2 ^ (^(^^t) 1/2^1/2^^^ obtain that \\G'/^{A^y^^\\ = || (^t) 1/2^1/2 Using this 
identity, the estimate above, we have for all and all u and all w = [G^^^'^u: 



ci > ||(At)^/^Gi/2||2 > n\-- J "-11 ^ ^^^^^^ ^ 



(^t)l/2gl/2^||2 ||(^t)V2^||2 



|(Gt)'/^ 



The estimate given above clearly implies that c^^{A^u, u) < (G^ii, u), and this is the lower bound 
in (5.2). The upper bound in (5.2) follows by interchanging the roles of G and A and basically 
repeating the same argument. □ 

The elementary results in the next proposition are used later in the proof of Lemma |5.5[ 



4 i _ 

Proposition 5.2. Let 9 G [0, 1] and define q{t; 6) = - — -(1 — - — -) and q{t; 9) = tq{t; 9). Then, 

9 -\- \ 6' + 1 
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(i) max q{t; 9) = 1; 

t<=ie,i] 

(ii) min q{t;e)= qi9; 9)= q{l; 9) ; 
te[6»,i] 

dq(l]9) 

(ill) — — — > (monotomcity). 
d9 

Proof. The proof of (i) and (ii) follow from the identity q{t;9) = 1 — {2t/{9 + 1) — l)^. The proof 
of (iii) is also straightforward and follows from the fact that 9 £ [0, 1] and hence 



de (^ + 1)2 v^ + i 

Next we derive estimates for the growth of the terms in a sequence, recursively defined using 
^(l; 9), which we use later to bound the convergence rate. 

4 t ~ 

Proposition 5.3. Let, 1 < c < 4 be a given constant, and q{t; 9) = ^(1 — -) and q{t; 9) = 

tq{t;9) (as in Proposition 5.2). Define, 

(5.3) 9i = l- 9k+i = -q0-]0k), for A; = 1,2,... 

c 

Then, the following are true for k = 1,2, . . .: 
(i) 4= - 1 ^ ^k+i <9k<l ; 



2 1 
(ii) 9k > max <j — — 1, 



/c 2A; - 1 + log A; _ 

Proof. The first item (i) follows from algebraic manipulations and the estimates given in Proposi- 
To show that 9k+i < 9^, we assume that 0^ > 2/-y/c— 1 (which is certainly true for k = 1. 



5.2 



tion 

To prove that ^fc+i > 2/^/c — 1 we observer that from 6*^ > 2/^/c — 1, the monotonicity property 
in Proposition 5.2 item (iii), implies that 

1 1 / 2 \ 2 

9k+i = -q{l; 9k)> -q{l;^-l] = ^-l. 

C C \ yJC J yJC 

Using again that ^fc > 2/ ^/c — 1 gives aso that 

The proof of the second item (ii) is a bit more involved. We prove this item by deriving an upper 
bound on (^fc = g^. Observe that, from the recurrence relation for 9^ we have 

(5.4) a+i = |(a + 2 + ^), ci = i- 

We first show that the faster growing sequence above is for c = 4. Indeed, let 

1 

Sk+i = Sfc + 2 H , si = 1. 

Sk 

A standard induction argument shows that 

Cfc < Sk, and 2k — 1 < Sk, V/c. 

Expand Sk by the recursive formula and we have 

fc-i ^ fe-i ^ 

Sfe = si + 2(A; - 1) + ^ - < 1 + 2(A: - 1) + ^ -— - < 2A: + In A: + 1, 
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which provides an upper bound of Cfc, and hence l/{2k + In A; — 1) is a lower bound of 9k- 

The following Lemma provides a spectral equivalence relation between G\. and B^^. 
Lemma 5.4. // Ai < X{B'^^Ak) < A2 and tq^it) > for Xi < t < X2, then 



(5.5) 



min{l, min tqk{t)} < ^ < max{l, max tqk{t)} 



Xl<t<X2 



yv:{v,l)=0, k = l,...,J -1. 



Proof. For any vector v, 

{qk{AkB,')v,B,'v) 
iAlv,v) 



{Qk{AlB^'Al){Al)^v,AlBj:'Al{Al)^v) _ {qk{Z)w,Zw) 



{aIv,v) 



{w, w) 
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where w = {A?)'^v and Z = A?B7^A?. Further, since Z has the same eigenvalues as BZ^A^, we 



conclude that 



min tqk(t) < 

\l<t<\2 



{qk{AkB7^)v,BZ^v) 



\l<t<X2 



This implies that for any x and y, 



x\ f {Y,l,Ak+,Yk+i)-' 



X 



a-'B,\k{AkB,')l \y 



'A fiYk+i^k+iYk+i. 



-1 



X 



{{Y^^^Ak+iYk+i)-^x, x) + a-\B^\{AkB^^)y, y) 
{{Y^^^Ak+iYk+i)-^x, x) + f7-i(A-iy, y) 

minjl, min maxjl, max tq{t)\ 

Xl<t<X2 Xi<t<\2 

and, hence, by using the definition of Gk and B^^, it follows that 



(5.6) 



minjl, min tGr(t)|, maxjl, max tq(t)} 

Xl<t<X2 Xi<t<X2 



□ 



Combining the above lemma with Theorem (4.11 ) the spectral equivalence between Bj^ ^ and A^., 
k = 1, . . . , J follows and is shown in the next Lemma. 

Lemma 5.5. Assume that the two level preconditioner Gk satisfies 

(5.7) {AkV,v) < {GkV,v) < Cg{AkV,v), \/v and k = 2, . . . , J. 
with constant Cg, such that 1 < < 4. Define 

(5.8) qk{t) =q{t,6k), 
where Ok are defined as 

01 = 1; e^^^ = ^q(l.ek) = -qkil). 



16 J. BRANNICK, Y. CHEN, J. KRAUS, AND L. ZIKATANOV 

Then, the following inequalities hold for all v : {v,l) =0 and k = 1, . . . , J. 



0k < \ \ ' < 1, 



(5.10) 



max 



1, 



1 



2k + lnk + l 



< 



Proof. We give a proof of (5.9) by induction. Clearly, for k = 1, B^^ = A\, and hence, (5.9) holds. 
We assume that the inequalities (5.9) hold for k = I and we aim to prove them for k = I + 1. For 
all V such that {v, 1) = we have 

{Bf_^\v,v) {Gl,v,v){Bf_^\v,v) 



Then, from (5.7), Proposition |5.1| and Proposition 5.4 (applied in that order) it follows that 
1 {Gj,-^v,v) {Bf_}^v,v) 

— < < 1, and minjl, min tqk(t)\ < < maxjl, max tqk(t)\. 

Next, by Proposition |5.2| and Proposition |5.3| we find that 



te[efe,i] 



1 {B^\v,v) 
= — min{l, min tqAt)} < < max{l, max tqi(t)} = 1. 



(5.11 



Finally, from the definition of i?^ ^ and ^ in terms of i?^ ^ and ^ , it immediately follows that 
{B,'v,v) _ {B-\Y,P)v,iY,P)v) 



< 



iAlv,v) 



{Y,P)v,{Y,P)v) 



<1, (t;,l) = 0. 



The proof of (5.10) follows from item (ii) in Proposition |5. 3 



□ 



The spectrum estimate (5.9) suggests that, Bj can be used as a preconditioner of a Conjugate 
Gradient method solving a linear system whose coefficient matrix is Aj. It also leads to the following 
convergence estimate of a power method. 

Theorem 5.6. Assume that there is a constant Cg such that 1 < Cg < 4 and {A^v, v) < {G^v, v) < 
Cg{AkV, v) for all v and k = 2, . . . , J . Then 



p((/ - nOd - BfA)) < mm ^j^^^^ 



< 1, 



where Hi is the £2 projection to the space of constant vectors. 



Proof. The proof is a directly application of the results in Lemma (5.5). 



□ 



A generalization of this estimate is given by assuming that Cg < for an integer m, in which 
case there exists an polynomial q{t) of order m — 1 such that a spectrally equivalent relation can 
be shown as 



m 



[m 



< -^^^^^ < 1, V-u : (-u,!) = and = 1,..., J, 



2-l)c„ (aIv,v 



which then implies that the power method preconditioned by the AMLI method using polynomial 
q{t) on all levels has a bounded convergence rate, as 

m^{cg - 1) 



p{{I-nr)iI-Bj'A)) < 



Cg{m? — 1) 
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For a matching on a hypercubic grid, as discussed above, the constant Cg approaches 4 asymptot- 
ically. Assume that the bound is given by Cg = 4, then a uniform convergence rate can not be 



proved by Theorem 5.6 since it requires that the two level spectrally equivalent constants on all 
levels must be less or equal to a common bound Cg which is strictly less than 4. This suggests us 
to find the best possible AMLI polynomials for the condition Cg = 4, and analyze how the AMLI 
convergence rate relates to the number of levels. 

Remark 5.7. An 1 — 1/logn type convergence rate can also be proven for the AMLI methods where 
the coarse partitioning consists of paths of m vertices where m> 2. 



6. Numerical results 

In the previous section, the convergence rate of two-level matching method was used to establish 
the convergence of the matching-based AMLI method. Here, a numerical implementation that is 
strictly a translation of this theoretical analysis is considered. Then, a simplified and more efficient 
variant of the method is developed and tested. 

To study the effectiveness of the algorithm and the sharpness of the theoretical estimates of 
its performance derived in the previous section, the method is applied as a preconditioner to the 
Conjugate Gradient iteration. In all tests, the stopping criteria for the PCG solver is set as a 10"^*^ 
reduction in the relative A norm of the error. The average convergence rate, r^, and the convergence 
rates computed by the condition number estimates obtained from the Lanczos algorithm and the 
AMLI polynomial, denoted by rg and r^, respectively, are reported. To reduce the effects of 
randomness in the numerical results, for each combination of testing parameters, the PCG method 
is run for five right hand sides computed by random left hand sides, and the convergence estimate 
that represents the worst case is reported. 



6.1. An exact implementation of the AMLI method. As a first test of the matching AMLI 
solver, it is applied to the graph Laplacian corresponding to 2- and 3-dimensional structured grids 
on convex and non-convex domains. The coarsening is obtained by applying matching only in a 
single direction on each level until the coarsest level is 1-dimensional, which is then solved using an 
LU factorization. The AMLI polynomial qk{t) on the A:-th level is determined by the theoretically 
estimated condition number, given by the recursive formula (5.4). The system YjfAf^Yk is solved 
exactly by an LU factorization on smaller grids or CG iteration down to 10~^ relative residual on 
larger grids of the hierarchy. 

Such AMLI method, which is designed to have all assumptions in Theorem 5.6 satisfied, is named 
"ordinary AMLI method." The results are reported in Table 6.1 and 6.2 and confirm that the actual 
convergence rate of the method, r^, and the condition number estimate, r^, match the theoretical 
estimate, that is, they both grow in accordance with the estimate = {^/k — l)/{-\/k + 1), where 
k grows logarithmically with respect to the grid size. 



6.2. Modified AMLI solver for matching. Next, a more practical variant of the matching 
AMLI preconditioner is developed. First, the exact A^Yk solvers are replaced by Richardson 
iterations with weights computed using the ii induced norm of these matrices, instead of the 
common choice of their largest eigenvalues. 

The lower order term In J = lnlog2 n in (5.4) is also dropped, since it is smaller than the term 2 J 
in (5.4) and is bounded by 4 for n = 2^*^. Another modification to the scheme is the choice of the 
scaling a in Lemma 4.4 away from 2. Numerical results suggest that a = 2 — 1/(2 log2 A'^), where N 
is the number of vertices on the graph, is usually a better scaling than the estimated bound a = 2 
used in the analysis. We use this choice for the structured mesh problems and for the unstructured 
problems the scaling is computed through a numerical method. 
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(a) Square domain with 
unknowns 



n 


k 


rk 


re 


ra 


128 


13.9 


0.58 


0.56 


0.54 


256 


16.0 


0.60 


0.59 


0.55 


512 


18.0 


0.62 


0.58 


0.57 


1024 


20.1 


0.64 


0.60 


0.60 


2048 


22.1 


0.65 


0.61 


0.61 



(b) L-shaped domain with 
(3/4)n^ unknowns 



n 


k 


rk 


re 


ra 


128 


13.9 


0.58 


0.56 


0.56 


256 


16.0 


0.60 


0.57 


0.59 


512 


18.0 


0.62 


0.57 


0.58 


1024 


20.1 


0.64 


0.59 


0.59 


2048 


22.1 


0.65 


0.60 


0.61 



Table 6.1. Results of the AMLI preconditioned CG method apphed to the 
graph Laplacians defined on 2D grids. 



(a) Cubic domain with 
unknowns 



n 


k 


rk 


re 


ra 


16 


16.0 


0.60 


0.55 


0.55 


32 


20.1 


0.64 


0.59 


0.59 


64 


24.2 


0.66 


0.62 


0.62 


128 


28.2 


0.68 


0.64 


0.64 



(b) Fichera domain with 
(7/8)n'^ unknowns 



n 


k 


rk 


re 


ra 


16 


16.0 


0.60 


0.55 


0.54 


32 


20.1 


0.64 


0.59 


0.59 


64 


24.2 


0.66 


0.62 


0.62 


128 


28.2 


0.68 


0.64 


0.64 



Table 6.2. Results of the ordinary AMLI preconditioned CG method applied 
to the graph Laplacians defined on 3D grids. 



In table 6.3 and 6.4 the convergence rate estimates of this approach applied to the same struc- 
tured problems are reported. Although some of the assumptions of the theory are violated by the 
method, its performance is similar to that of the approach considered in the previous tests. 

Remark 6.1. A more practical strategy is to use a numerical method, e.g., a Lanczos algorithm 
with an AMLI preconditioner on the k-th level, to estimate the smallest eigenvalue of Bj^^Aj^, which 
is then used to determine the AMLI polynomial on the k + 1-th level. Numerical tests show that 
such strategy results faster convergent AMLI methods than that defined through recursive formula 



(5.4), at a cost of more complicated setup phase. This strategy usually provide a significant speed 



up for 3- or higher dimensional structured problems. 

6.3. On unstructured grids. Finally, tests of this AMLI preconditioned Conjugate gradient 
method applied to the graph Laplacian defined on more general graphs, coming from unstructured 
meshes resulting from triangulations of a 2-dimensional grid on a square domain, or a 3-dimensional 
grid on a cubic domain, are considered. The unstructured grid is generated by perturbing grid points 
of a structured grid by a random vector of length h/2, where h is the mesh size of the original struc- 
tured grid, followed by a Delaunay triangulation. Then, a random matching is applied recursively 
to generate a multilevel hierarchy with (log2 N)/2 levels. The 3-dimensional unstructured grids are 
generated in a similar way and the multilevel hierarchy is constructed accordingly by the random 
matching algorithm. 



The results of these tests are reported in Table 6.5 and 6.6 For the results on the left of 
these tables, the Y^AkYk block of the two-level preconditioner is solved to high accuracy, which is 
practical since this operator is proven well conditioned even for unstructured grids. The recursive 
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(a) Square domain with 
unknowns 



n 


k 


rk 


re 


ra 


128 


13.0 


0.57 


0.59 


0.54 


256 


15.0 


0.59 


0.62 


0.58 


512 


17.0 


0.61 


0.64 


0.59 


1024 


19.0 


0.63 


0.65 


0.63 


2048 


21.0 


0.64 


0.65 


0.65 



(b) L-shaped domain with 
(3/4)n^ unknowns 



n 


k 


rk 


re 


ra 


128 


13.0 


0.57 


0.58 


0.56 


256 


15.0 


0.59 


0.60 


0.56 


512 


17.0 


0.61 


0.62 


0.57 


1024 


19.0 


0.63 


0.64 


0.62 


2048 


21.0 


0.64 


0.69 


0.67 



Table 6.3. Results of the modified AMLI preconditioned CG method apphed 
to the graph Laplacians defined on 2D grids. 



(a) Cubic domain with (b) Fichera domain with 

unknowns {7/8)n'^ unknowns 



n 


k 


rk 


re 


ra 




n 


k 


rk 


re 


ra 


16 


15.0 


0.59 


0.50 


0.42 




16 


15.0 


0.59 


0.49 


0.49 


32 


19.0 


0.63 


0.54 


0.49 




32 


19.0 


0.63 


0.54 


0.50 


64 


23.0 


0.65 


0.57 


0.52 




64 


23.0 


0.65 


0.57 


0.56 


128 


27.0 


0.68 


0.59 


0.56 




128 


27.0 


0.68 


0.57 


0.60 



Table 6.4. Results of the modified AMLI preconditioned CG method applied 
to the graph Laplacians defined on 3D grids. 



formula (5.8) is used to derive the polynomials used in the AMLI cycles, and the scaling constants 



are computed using 



{P^AkPkh 
ak = max — f- — - 

ij^j {Ak+l)ij 



which ensures that the upper bound in (5.9) is always 1, which in turn guarantees that the AMLI 



method, as a preconditioner for the CG method, is always positive semi-definite. Because that the 



AMLI polynomials, constructed according to (5.8), is negative when t > 1. Assume that the scaling 



constant ak is smaller than the value suggested above, then there exists a v such that 

{Glv,v)>{Alv,v), 
it. 



which makes it possible that {Bf^^v, v) > {Af,v, v). Assume that happens, the matrix B^^q{AkBjl^^) 
becomes indefinite which in turn makes -Bfc +i indefinite. 

For the results on the right of Table 6.5 and 6.6, the solve of the AkYk block is replaced by 



one Richardson iteration, and the AMLI polynomials are constructed based on (5.8) without the 



lower order term In k. The asymptotic convergence rates are again close to the expected convergence 
rates obtained from the AMLI polynomials Further, the actual convergence rates are usually better, 
especially for the method that uses more accurate solves for the VjFAkYk blocks, as opposed to the 
one that uses a single Richardson iteration. 
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(a) Ordinary AMLI (b) Modified AMLI 



n 


k 






ra 




n 


k 


rk 






128 


16.0 


0.60 


0.70 


0.58 




128 


15.0 


0.59 


0.70 


0.70 


256 


18.0 


0.62 


0.72 


0.54 




256 


17.0 


0.61 


0.71 


0.70 


512 


20.1 


0.64 


0.74 


0.63 




512 


19.0 


0.63 


0.72 


0.72 


1024 


22.1 


0.65 


0.75 


0.65 




1024 


21.0 


0.64 


0.73 


0.73 


2048 


24.2 


0.66 


0.76 


0.67 




2048 


23.0 


0.65 


0.75 


0.75 



Table 6.5. Results of the CG method preconditioned by variants of the matching 
AMLI methods applied to the graph Laplacian defined on 2D unstructured grids of 
size n^. 



(a) Ordinary AMLI 



(b) Modified AMLI 



n 


k 




re 




16 


18.0 


0.62 


0.65 


0.48 


32 


22.1 


0.65 


0.67 


0.55 


64 


26.2 


0.67 


0.70 


0.62 


128 


30.3 


0.69 


0.74 


0.60 



n 


k 


Tk 






16 


17.0 


0.61 


0.59 


0.55 


32 


21.0 


0.64 


0.63 


0.58 


64 


25.0 


0.67 


0.65 


0.62 


128 


29.0 


0.69 


0.67 


0.65 



Table 6.6. Results of the CG method preconditioned by variants of the matching 
AMLI methods applied to the graph Laplacian defined on 3D unstructured grids of 



size n'^ 



7. Conclusions 

An algebraic formula for estimating the convergence rate of an aggregation-based two level 

method is derived, and it is shown that the formula can be used to obtain sharp estimates of 
the convergence rates in the special case where matching is used. With the use of geometric 
information, a sharp bound of the two-level method is derived. The nearly optimal convergence 
and complexity of the multilevel method that uses AMLI cycles is also established. The reported 
numerical tests illustrate the sharpness of the theoretical estimates. Moreover all the theoretical 
results can be generalized to aggregates of general size and, hence, can be used to study an approach 
which combines aggressive aggregation with AMLI cycles, which should result in a fast and memory 
efficient solver for graph Laplacians. Development and analysis of such a scheme and one that uses 
more general smoothers are subject of on-going research. 
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